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A detailed derivation of the recently proposed time-dependent numerical renormalization-group 
(TD-NRG) approach to nonequilibrium dynamics in quantum impurity systems is presented. We 
demonstrate that the method is suitable for fermionic as well as bosonic baths. A comparison with 
exact analytical results for the charge relaxation in the resonant-level model and for dephasing in 
the spin-boson model establishes the accuracy of the method. The real-time dynamics of a single 
spin coupled to both types of baths is investigated. We use the TD-NRG to calculate the spin 
relaxation and spin precession of a single Kondo impurity. The short- and long-time dynamics 
is studied as a function of temperature in the ferromagnetic and antiferromagnetic regimes. The 
short-time dynamics agrees very well with analytical results obtained at second order in the exchange 
coupling J. In the ferromagnetic regime, the long-time spin decay is described by the scaling variable 
x — 2pF.J(T)Tt. In the antiferromagnetic regime it is governed for T < Tk by the Kondo time 
scale 1/Tk ■ Here pF is the conduction-electron density of states and Tk is the Kondo temperature. 
Results for spin precession are obtained by rotating the external magnetic field from the x axis to 
the z axis. 

PACS numbers: 03.65.Yz, 73.21. La, 73.63.Kv, 76.20.+q 
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I. INTRODUCTION 



The decoherence and relaxation of an impurity spin is 
a classic problem in condensed mater physics. More then 
30 years ago, Langreth and Wilkins developed a theory 
of spin resonance in dilute magnetic alloys. 1 The princi- 
pal objective was to derive Bloch-like equations for the 
paramagnetic resonance using the Kadanoff-Baym^ and 
Keldysb-. techniques to nonequilibrium. This included 
the description of spin precession and spin relaxation of 
a finite concentration of magnetic impurities weakly cou- 
pled to a metallic host. The recent advent of quantum- 
dot devices and their potential application to quantum 
computing has renewed interest in the spin dynamics of a 
single impurity. In contrast to real magnetic impurities, 
quantum dots can be controlled in exquisite detail, and 
can be tuned at will from weak coupling to the Kondo 
regime. Perturbative approaches, such as that of Lan- 
greth and Wilkinsji fail to describe the strong correlation 
that develop in latter regime. This calls for new theoret- 
ical approaches to nonequilibrium dynamics, suitable for 
treating the strongly correlated state. 

In this work, we focus on the dynamics of a single spin 
coupled to an infinitely large environment of noninteract- 
ing particles, in order to investigate spin-relaxation phe- 
nomena at all temperatures and coupling regimes. Re- 
laxation of a single spin is the simplest example for real- 
time dynamics in a quantum impurity system. Other 
examples might be qubits, coupled quantum dots, or bi- 
ological donor-acceptor molecules. In such systems one 
is generically interested in the dynamics of a finite sub- 
system interacting with an infinitely large environment. 
The nature of the environment may vary from one real- 
ization to another. It can consist of a fermionic bath, 



as in single-electron and single-molecule transistors, or a 
bosonic bath, as in the spin-boson model. In certain cases 
a combined fermionic-bosonic bath might be in order. 

The immense difficulty of treating the real-time dy- 
namics of quantum impurity systems stems from the need 
to track the full time evolution of the density operator 
of the entire system — environment plus impurity. The 
Kadanoff-Baym and Keldysh techniques^ provide an el- 
egant platform for perturbative expansions of the density 
operator. In general, however, perturbative approaches 
are plagued by the infra-red divergences caused by de- 
generacies on the impurity, making them inadequate for 
tackling quantum impurities. Here we take an alternative 
approach to the real-time dynamics based on Wilson's 
numerical renormalization-group (NRG) method^ 

The NRG is a very powerful tool for accurately cal- 
culating equilibrium properties of arbitrarily complex 
quantum impurities. Originally developed for treating 
the single-channel Kondo Hamiltonianj* this nonpertur- 
bative approach was successfully extended to the An- 
derson impurity model (symmetries and asymmetrio&) , 
the two-channel Kondo Hamiltonian, 7 - 8 different two- 
impurity clusters (SiiSiiiiiSiiSiii and a host of related zero- 
dimensional problems. Recently, we extended the ap- 
proach to a certain class of time-dependent problems 
where a sudden perturbation is applied to the impu- 
rity at time t = 0M> Similar to the equilibrium NRG, 
the time-dependent NRG (TD-NRG) can be applied to 
all coupling strengths and is not confined to the weak- 
coupling regime. It is capable of accessing all time scales, 
including arbitrary long as well as arbitrary short times. 
These appealing properties make the TD-NRG a power- 
ful new approach for studying nonequilibrium dynamics 
in quantum-impurity systems. 

In this work, we present the complete details of the TD- 



NRG, and apply it to the spin dynamics of the Kondo 
model. A comprehensive analysis is presented for the 
spin relaxation and spin precession as a function of tem- 
perature, magnetic field and coupling strength. Both fer- 
romagnetic and antifcrromagnetic couplings are consid- 
ered. In the Kondo regime, spin dynamics is governed by 
the Kondo time scale tx = 1/Tk, where Tk is the Kondo 
temperature of the system. In the ferromagnetic regime, 
the underlying time scale is strongly dependent on tem- 
perature, reflecting the fact that the effective equilibrium 
coupling flows to zero with decreasing temperature. 



A. Preliminaries 

The idea of applying the NRG to nonequilibrium dy- 
namics dates back to the work of Costi^ In the spirit 
of the equilibrium NRG, Costi evaluated the nonequilib- 
rium spectral functions using an individual Wilson shell 
for each frequency interval. However, as already recog- 
nized by Costifi^ a full multiple-shell evaluation is ulti- 
mately required for the correct description of nonequi- 
librium dynamics. The physical reason is simple: When 
expanded in terms of the eigenstates of the perturbed 
Hamiltonian, the initial state of the unperturbed system 
contains contributions from all energy scales. This calls 
for an adequate coupling of low- and high-energy scales 
absent in the equilibrium NRG. 

To circumvent this problem, we construct a complete 
basis set of the Wilson chain using the NRG eigen- 
states. Implementing a suitable resummation procedure, 
we track all states that contribute to the time evolution 
of the observables of interest. 15 Explicitly, we have shown 
that the time evolution of a general local operator O can 
be written as 



N trun 



(1) 



where Ey™ an( i ^T are the dimension-full NRG eigen- 
energies of the perturbed Hamiltonian at iteration m < 
N, OJ™ s is the matrix representation of O at that iter- 
ation, and p r s e ^(m) is the reduced density matrbsiS in 
the basis of perturbed Hamiltonian. The restrictive sum 
over r and s requires that at least one of these states is 
discarded at iteration m. In context of the equilibrium 
NRG, the reduced density matrix has been already used 
for the calculation of spectral functions^ 

Equation (JIJ is the centerpiece of the TD-NRG ap- 
proach and will be heavily using throughout this paper. 
Below we present a detailed discussion of its derivation 
and implementation. Here we only wish to emphasize 
the following points, (i) The reduced density matrix oc- 
curs naturally in our formulation, (ii) It does not fol- 
low a unitary time evolution and, therefore, contains in- 
formation on dissipation and decoherence. (iii) Equa- 
tion Q arises from summation over the complete many- 
body Fock space of the Wilson chain. No truncation of 



states is involved. 



B. Plan of the Paper 

As implied above, this paper has two main objectives: 
an indepth presentation of the TD-NRG and a compre- 
hensive analysis of the spin relaxation for a single spin 
coupled to a conduction band. Accordingly, the remain- 
der of the paper divides into two major parts. 

The following three sections are devoted to an exten- 
sive exposition of the TD-NRG. This includes a detailed 
derivation of the approach, along with benchmark results 
for both a fermionic and a bosonic bath. Section III Bl 
contains a general formulation of TD-NRG for any ini- 
tial density operator po of the system. In particular, 
Eq. Q is derived and an explicit definition of the re- 
ducted density matrix // ed (m) is given. In Sec. Ill Gl we 
discuss how the reduced density matrix is calculated in 
practice in the generic case where po is the equilibrium 
density operator for the unperturbed system. The com- 
plete TD-NRG algorithm is outlined in turn in Sec. Ill Dl 
In sections IlIII and llVl we compare the TD-NRG to ex- 
act results for two distinct models: the resonant-level 
model and a restricted version of the spin-boson model. 
These two test cases allow us to establish the accuracy of 
the TD-NRG on all time-scales, and to investigate how 
well a discretized finite-size system can be used to mimic 
nonequilibrium dynamics in a continuous bath. 

Following this detailed exposition of the TD-NRG, we 
proceed in sections IVI and IVT1 to the spin dynamics in the 
Kondo model. In Sec. [V] we introduce the model and 
present some basic analytical considerations. These in- 
clude an explicit calculation of the short-time response 
up to second order in the dimensionless exchange cou- 
pling. Numerical results for arbitrarily long time scale 
are covered in turn in Sec. I VII both for the ferromagnetic 
and antiferromagnetic regimes. Different configurations 
of the magnetic field are also considered. Through com- 
parison with the analytical results of Sec. E] we are able 
to establish the accuracy of the TD-NRG results at short 
time scales. The true power of the numerical solution 
lies, however, in the long-time behavior not accessible by 
perturbative techniques. 



II. THEORETICAL FORMULATION 

A. Goal of the method 

In the present section we formulate the TD-NRG ap- 
proach in general terms, before turning to concrete ex- 
amples in sections IIIII IIVI and IVII Hereafter we assume 
that the initial density operator at time t = is known, 
when a sudden perturbation A7i is switch on. For t > 0, 
the system evolves in time with the full Hamiltonian 
a? = H l + AH. where TV 1 is the Hamiltonian of the 



unperturbed system at time t < 0. Denoting the ini- 
tial density operator by p , the expectation value of any 
operator O at time t > is given by 



0(t) = Tr{p(t)0} 



(2) 



where p(t) is the corresponding time-dependent density 
operator. 

Our first goal is to show that expectation value 0(t) of 
any local operator O can be rewritten in terms of a sum 
over the contributions of a sequence of reduced density 
matrices p red (t)m^H The reduced density matrices are 
generated by systematically tracing out all environment 
degrees of freedom. Hence the unitary time evolution of 
the density matrix of the system as a whole is lost, giv- 
ing rise to dissipation and decoherence. Although our 
derivation makes explicit reference to the NRG proce- 
dure, it might be applicable to other methods as well. 

Two key ingredients underlie our approach: (i) The 
identification of a complete basis set of the many-body 
Fock space of the Wilson chain based on the NRG eigen- 
states at the different iterations; (ii) Expectation values 
are obtained by explicitly tracing over this complete basis 
set. This marks a significant departure from the conven- 
tional renormalization-group concept, where the relevant 
physical information is contained in the kept states. Here, 
the discarded states are solemnly responsible for the time 
evolution. This also solves the problem of overcounting 
excitations encountered by Costii^S, Each excitation con- 
tributes only once to expectation values, at that iteration 
where the corresponding state is discarded. For the case 
where po either projects onto a single state or represents 
the full equilibrium density operator of the unperturbed 
Hamiltonian, we are able give a closed analytical expres- 
sion for the reduced density matrices. 



Real-time evolution in quantum impurity 
systems 




m,N 



FIG. 1: The full Wilson chain of length N is divided into 
a sub-chain of length m and the "environment" R m ,N- The 
Hamiltonian TL m can be viewed either as acting only on the 
sub-chain of length m, or as acting on the full chain of length 
N, but with the hopping matrix elements t m , ■ ■ ■ , ijv_i all set 
to zero. The former picture is the traditional one. In this 
paper we adopt the latter point of view. 



A > 1 A The continuum limit is recovered for A — > 1 . Us- 
ing an appropriate unitary transformation^ the Hamil- 
tonian is mapped onto a semi-infinite chain, with the 
impurity coupled to the open end. The Nth. link along 
the chain represents an exponentially decreasing energy 
scale: Dm ~ A~ N / 2 for a fermionic batbJ and Dn ~ A~ N 
for a bosonic bathi 9 Using this hierarchy of scales, the 
sequence of finite-size Hamiltonians Hn for the A-site 
chain^i is solved iteratively, discarding the high-energy 
states at the conclusion of each step to maintain a man- 
ageable number of states. The reduced basis set o{Hn so 
obtained is expected to faithfully describe the spectrum 
of the full Hamiltonian on a scale of Dn, corresponding 
to the temperature Tn ~ Dn~ 

Due to the exponential form of the Boltzmann factors, 
the reduced NRG basis set of TLm is sufficient for an ac- 
curate calculation of thermodynamic quantities at tem- 
perature Tpj. This is no longer the case once the system 
is driven out of equilibrium. Since the nonequilibrium 
dynamics involves all energy scales exceeding T/v, and 
in the absence of a general criterion as to which states 
contribute to the dynamics, a complete basis set of the 
Fock space Tn of Hn is required. In the following we 
identify such a complete basis set of Tn, composed of 
approximate eigenstates of Hn- 



The Hamiltonian of a quantum impurity system is gen- 
erally given by 



H = 7Yb 



ath 



■Hi, 



Ht, 



(3) 



where 7ibath models the continuous bath, 7ii mp represents 
the decoupled impurity, and 7i m i X describes the coupling 
between the two subsystems. The entire system is char- 
acterized at time t = by an arbitrary density matrix 
po, when a sudden, time-independent perturbation AH 
is switched on: H(t > 0) = H l + AH = H f . For t > 0, 
the density operator evolves according to 



m 



e-™'%e i7lft 



(4) 



In equilibrium, one can solve such a quantum impurity 
system very accurately using the NRG. At the heart of 
this approach is a logarithmic discretization of the con- 
tinuous bath, controlled by the discretization parameter 



1. Complete Basis Set 

There are two possible ways to envision the iterative 
NRG solution of the iV-site chain. In the traditional pic- 
ture one starts from a core cluster that consists of the im- 
purity degrees of freedom and the n — site, and enlarges 
the chain by one site at each NRG step. Alternatively, 
one can view the NRG procedure as starting from the full 
chain of length N, but with the hopping matrix elements 
set to zero along the chain. At each successive step an- 
other hopping matrix element is switched on, until the 
full Hamiltonian Hn is recovered. In this latter picture, 
to be adopted below, the entire sequence of Hamiltoni- 
ans H m with m < N act on the Fock space of the A-site 
chain. Accordingly, each NRG eigen-energy of H m has 
an extra degeneracy of dS m ', where d is the number 
of distinct states at each site along the chain. The extra 
degeneracy stems from the N — m "environment" sites 



at the end of the chain, depicted in Fig. ^ 

Let us elaborate on the formal connection between 
the two pictures presented above. Let {|r;m)} label the 
NRG eigenstates of Ti m when acting on the TO-site chain, 
and let E™ denote their corresponding eigen-energies. 
Enumerating the different configurations of site i by 
{a;}, each of the tensor-product states |r, a m +i, • • • , a^) 
with arbitrary a m +i, • • • , ajv is then a degenerate eigen- 
state of H m with energy E™, when acting on the full 
iV-site chain. For brevity, we introduce the short- 
hand notation \r,e;m), where the "environment" vari- 
able e = {a m +i, • • • , chat} encodes the N — m site labels 
a m+ i, • • • , a at. The index to is used in this notation to 
record where the chain is partitioned into a "subsystem" 
and an "environment" (see Fig. Q. 

In the conventional NRG picture, the eigenstates of 
TL m +i are constructed from the tensor-product states 
{|r, a m +i; 771)}, corresponding to enlarging the chain by 
one extra site. A unitary transformation U relates 
the new eigenstates \r';m + 1) of H m +i to the basis 
{\r,a m+1 ;m)}: 

\r';m+l) = ^ U r',ra m +i \r, a m +i; to) . (5) 

r,a m + 1 

An alternative notation is given by 

|r';m + l) = ^ P"> l A a m+i]\r,a m+ i;m) , (6) 

r,a m + 1 

where the matrix elements P r * ^ r [a m +\\ are identified with 
the corresponding matrix elements of the unitary trans- 
formation U: P™ r [a m +i] = U r ' tram+1 . Successively 
applying the recursion relation of Eq. @, the eigen- 
states of H.n can be formally viewed as matrix-product 
states^ generated by consecutive application of the ma- 
trices P™ r {a m+1 }. 

Note that the transition itself from the tensor-product 
states {\r,a m +i;m}} to the eigenstates {\r';m + 1)} in- 
volves no truncation of the Fock space. However, since 
the dimension of the Fock space grows as d N , a complete 
basis set is not manageable for any practical chain length 
of order N ~ 100. In the equilibrium NRG, high-energy 
states are thus discarded after each iteration, as these do 
not contribute to the equilibrium density matrix. This 
latter statement is guarantied by the hierarchy of scales 
along the Wilson chain. It would not apply to any or- 
dinary tight-binding chain with constant hopping matrix 
elements. 

Consider now the first iteration to at which states are 
discarded. In order to keep track of the complete basis 
set of the TV-site chain, we formally divide the eigenstates 
|r, e;ro) into two distinct subsets: the discarded high- 
energy states {\l , e; m) dts} and the retained low-energy 
states {\k , e; m) k P } ■ In the course of the NRG, only the 
latter states are used to span the Hamiltonian 7i m +i 
within the reduced subspace {\k, a m +i,e';m)}. Note, 
however, that the sum of both subsets still constitutes 
a complete basis set for the JV-site chain. Repeating this 



procedure at each subsequent iteration, we recursively 
divide the retained subset into a discarded part and a 
retained part. Then, the collection of all discarded eigen- 
states \l,e;m)dis together with the eigenstates of the fi- 
nal NRG iteration N combine to form a complete basis 
set for the entire Fock space Tn- Regarding all eigen- 
states of the final NRG iteration as discarded, one can 
formally write the Fock space of the iV-site chain in the 
form T N = span{|Z, e; m) dis }. 

We stress that all states are retained in the course of 
this construction. Not a single state is eliminated. Since 
{|Z, e; m)dis\ constitutes a complete basis set of J-n, then 
the following completeness relation obviously holds: 



E 



2_j \ l i e > m )dis dis(l,e;'i 



= 1 



(7) 



l.e 



Here the summation over to starts from the first iteration 
TOmin at which a basis-set reduction is imposed. Typical 
values of TO m i n are 4-5 for a spin-degenerate conduction 
bath. The summation indices I and e implicitly depend 
on to. The identification of this complete basis set, natu- 
rally generated by the NRG, is one of the two key ingre- 
dients of our method. All traces will be carried out below 
with respect to this basis set. Hence, the evaluation of 
time- dependent expectation values involve no truncation 
error. Note, however, that we made no reference to a 
particular Hamiltonian TL in constructing the basis set 
{\l,e;m)dis}- Below we shall make use of two distinct 
basis sets of this form, one for the Hamiltonian TL? and 
another for the unperturbed Hamiltonian TL 1 . 

Another useful identity to be used below pertains to 
the subspace retained at iteration to, {\k,e;m)kp\- To 
this end, we note that the sum in Eq. can always be 
divided into two complementary parts 1~ and 1+ : 

m 

l m = E ^2\l',e';m')dis dis(l',e';m'\ , (8) 

TV 

X m = E $> , .e / ;ro'>*.«K.<l , ,e';m , |. (9) 

m'— m+1 I' ,e' 

(For to = N only 1~ exists.) The completeness relation 
therefore becomes 



1 



-^m, i -*-r 



(10) 



What do the operators 1~ and 1+ represent? The 
operator 1~ projects onto the subspace of all discarded 
states up to and including iteration to. The operator 
1^ projects onto the complementary subspace retained 
at iteration m. An alternative way of writing the latter 
projection operator is in terms of the states retained at 
iteration to, {\k,e;m)kp\- Explicitly, 1+ is given by 



£ - E 

k,e 



k, 



) kp kp 



(fc, e; m\ 



(ii) 



Equations (|5J and (|llfl provide us with an important 
connection between summations over retained and dis- 
carded states. These equivalent representations of 1+ 
will be frequently used later on. 



2. Time-evolution of local expectation values 

We are now in position to formally evaluate the time- 
dependent expectation value of any local operator O at 
time t > 0. As indicated in Eq. Q, we need to explicitly 
carry out the trace over p(t)0. This is most conveniently 
done using the basis set {\l , e; m) dis} introduced in the 
previous section: 



JV 



°W = E ^ ^a,e;7Ti|p(t)0|Z,e; 



m)di S 



(12) 



771 — m m i n l.e 



Inserting Eq. (|1L)|) in between the operators p(t) and O 
in Eq. (U2J) yields 



N 



°W= E I]^,e;m|p(*)(l-+l+)6|/,e; 



TO><ji 



m— rri m i n /,, 



(13) 
Using Eqs. (JSJ) and (fTTfl for l m and 1+, respectively, we 
obtain 6{t) = 0_(t) + O (t) + 0+{t), where 

JV m— 1 

°-(*) = E E EE^f 1 ' 6 '™'!^-^)* 

m>m m i n m'— m m i n Z,T e,e f 

x d j s (Z, e; m|/5(?;)|/', e'; m')^ (14) 



and 



JV 



°o(t) = J^ EE di s (^e;m|p(t)|/',e';m)di s 

m=m m i n Z,P e,e' 

Xdis(l\ e';m\0\l,e;m)dis (15) 

stem for the projection operator 1~ , and 

JV-l 

°+(*) = E ^2^2 dis(l,e;m\p(t)\k,e';m) kp 



m=m TaiTl l y k e,e' 

x kp (k, e';m\d\l,e;m)di 



(16) 



originates from 1+. 

Whereas the terms Oo(t) and 0+(i) involve the sum- 
mation over a single iteration variable to, 0_ (t) contains 
two such summations: one over to and another over m' . 
Rearranging the two sums according to 

N 777-1 7V-1 IV 

E E - E E . d7) 

m=m m i n +l m/— m m i n tti 7 — 7n m i n tti— m' + l 

utilizing the identity 

TV 

E Ei'' e;m ^ isdis ^' e;m i = 1 ™' ' ( 18 ) 

777— m' + l /.e 



and replacing 1+, by the representation of Eq. (|11|) we 
obtain 

N-l 

°-W = E SZ)*p(*.e;m / |^(*)|f / ,e / ;m / )d«i 

rn 7 — m m i n fe,i' e,e' 

x d ,- a (f', e'; m'|6|fc, e; m')fc P . (19) 

Combining Eqs. ifTB - ]!. (fTf))) and fli^l then yields 

TV trun 

°(*) = E EE ^H^Ke',™) 



x (r, e'; m|0|s, e; to) 



(20) 



Here the restricted sum ^Z r ™ ra implies that at least one 
of the states r and s is discarded at iteration to. Terms 
where both states are retained contribute to the sum at 
some later iteration ml > to. 

So far we made no assumption about the operator O, 
nor have we committed ourselves to a particular Hamilto- 
nian in writing the basis set {\l, e;m)dis}- In the follow- 
ing we restrict ourselves to local operators O. By local 
we mean that the operator acts on degrees of freedom 
that reside either on the impurity or on close-by sites to 
for which all states are still available (i.e., m < m m i n ). 
This restriction is rather weak, since all operators in the 
vicinity of the impurity fulfill this requirement. The ma- 
trix elements of a local operator O are diagonal in and 
independent of the environment degrees of freedom: 

(r,e;m\d\s,e';m) = 6 e ,e'0™ s , (21) 

This allows us to write Eq. (|20|l in the form 

N trun 

°W= E EE ™^ s ' e ' m lpWk,e;TO). (22) 

m— m m i n r,s e 

As for the states {\l,e;m)dis}, unless stated otherwise 
we work hereafter in the NRG basis set generated for the 
perturbed Hamiltonian Tif . This choice is motivated by 
the time dependence of pit), which is governed by W . 
Indeed, using Eq. J3J and resorting to the standard NRG 
approximation— H N \k, e; to) ~ E™\k, e; m) one has 



; m\p(t)\r, e; to) ss e 



i(E™-E™)t 



s,e;m\p \r,e;m) 



(23) 

Upon inserting Eq. I|23(l into Eq. 122(1 . the environment 
variable e enters only through the matrix element of po- 
Introducing the reduced density matrisii^ 



PT,r( m ) = E( s ' e; m l/Sok, e; 



(24) 



one arrives at the final result for the time evolution of 
0(t) at t > 0^ 

N truri 

°W = E E e l(E ™- ET)t OZ P\ cd r {m) . (25) 



Several comments should be made about Eq. (|25|l . 
First, no assumption was made about the initial density 
operator p . It can either project onto a particular state 
or represent the equilibrium density operator for TC . It 
can even stand for a density operator that has already 
evolved in time subject to some previous dynamics. Sec- 
ond, all states of the finite-size Fock space are retained in 
Eq. I|25[l and all energy scales are explicitly taken into ac- 
count. No basis set reduction is imposed. This should be 
contrasted with recent time-dependent extensions of the 
density-matrix renormalization group i2i*22iSi2li2^ In the 
so-called TD-DMRG, the ground state is evolved in time 
using a Trotter decomposition of the time-evolution op- 
erator. This introduces an accumulated error that grows 
linearly in time. Here 0(t) is evaluated independently for 
each time t > 0, and is thus free of accumulated errors. 

This does not mean that the Eq. (|25f) is error- free. 
Two approximations enter the TD-NRG: (i) a discretized 
finite-size representation of the continuous bath, and 
(ii) the conventional NRG approximation H N \k : e\m) rs 
E™\k, e; m). Indeed, the latter approximation is the only 
one invoked in describing the time evolution of expecta- 
tion values on the iV-site chain. In particular, Eq. 12">l) is 
exact for t — ► + . As shown by Wilson^ the associated 
error in thermodynamic quantities is perturbative and 
small, a consequence of the separation of scales along the 
Wilson chain. To assess the accuracy of this approxima- 
tion in the present context we have compared Eq. I|25|) 
to the exact time evolution on sufficiently short chains 
where all eigenstates of H? can be computed and stored 
in memory. To facilitate a meaningful comparison, only 
a deliberately small number of the states were kept in 
the course of the NRG iterations. A surprisingly good 
agreement was found between the exact results and the 
TD-NRG up to very long times, establishing the accu- 
racy of the accepted NRG approximation in the present 
context. 



A more significant source of error is due to the dis- 
cretized finite-size representation of the continuous bath. 
There are two aspects to this approximation. Due to the 
limited energy resolution at low energies, Eq. I|25(l gener- 
ally looses accuracy for t 3> 1/Djv- Since the chain length 
N is chosen such that Dn ~ T, one expects then a lose 
of accuracy for t ^> 1/T. It should be noted, however, 
that the NRG can access arbitrarily low temperatures, 
implying that arbitrarily long time scales can be reached 
for T — > 0. We demonstrate this important point later 
on. At the same time, a continuous spectrum is expected 
to be vital for relaxation to the exact thermodynamic ex- 
pectation value with respect to Ji' . This constitutes a 
fundamental obstacle for any solution, however accurate, 
of a finite-size system. As detailed in subsection III El we 
can largely overcome this obstacle by (i) averaging over 
different realizations of the Wilson chain using Oliveira's 
z-trickjSi and (ii) resorting to a Lorentzian broadening 
of the NRG levels. 



C. Calculation of the reduced density matrix 

So far, our discussion applied to a general local oper- 
ator O and to an arbitrary initial density operator p . 
The time evolution of 0(t) was shown to be determined 
by a sequence of reduced density matrices p rcd , in which 
the environment degrees of freedom are traced out it- 
eratively. While formally applicable to any /3o, practical 
calculations may prove more restrictive. Here we address 
the calculation of the reduced density matrix. 

Typically, po has a simple representation with respect 
to the eigenstates of the unperturbed Hamiltonian Ti % . 
For example, starting from thermal equilibrium po is 
equal to e-F^/Zi, where Z, L = Tr{e- fm '} is the un- 
perturbed partition function. Thus, the reduced density 
matrix is most conveniently expressed with respect to 
the NRG eigenstates of the unperturbed Hamiltonian, 
i.e., when the states {\l,e;m}} in Eq. l(^|) relate to Ti l . 
Note, however, that the reduced density matrix which 
enters Eq. (|2*3|l was explicitly constructed for the NRG 
eigenstates of the full Hamiltonian Ti^ . In general, there 
is no simple relation between the two reduced density 
matrices, as they involve traces over different environ- 
ments. Still, one can convert between the two entities in 
the generic case where the system is evolved in time from 
thermal equilibrium. Below we explain in detail how the 
reduced density matrix of Eq. (|24[1 is calculated for this 
generic case. 



1. General considerations 

Let us assume for the time being that the reduced den- 
sity matrix is know with respect to the NRG eigenstates 
of the unperturbed Hamiltonian TC 1 . Namely, when the 
states {\l,e;m}} that enter Eq. (|24|l correspond to TC" . 
Our goal is to compute the reduced density matrix with 
respect to the NRG eigenstates of the full Hamiltonian 
H.f . Technically, this requires working with two distinct 
sets of NRG states — one for Ti 1 and the other for H? . To 
simplify the notation as much as possible, we distinguish 
the two sets of states by their labels. Throughout this 
subsection we designate the NRG states pertaining to the 
unperturbed Hamiltonian 7i l by indices carrying the sub- 
script i. For example, the state \li,ei\m). NRG states 
corresponding to TV will be labeled, as before, with plain 
indices having no additional subscripts as in \l, e; to). We 
stress that ei and e label the same environment degrees of 
freedom, despite the difference in appearance. Similarly, 
we reserve the notation p r s ed (m) for the reduced density 
matrix as defined in Eq. I)24|) with respect to the states 
of the full Hamiltonian Ti^ . The reduced density matrix 
with respect to the states of the unperturbed Hamilto- 
nian W is denoted (and defined) by 



„rcd,0 



(m) 



{s i ,e i \m\p \r l ,e l ]m) 



(26) 



As shown in the previous subsection, one can write the 



completeness relation 1 = 1~ + 1+ for any Hamiltonian 
H. The projection operators 1~ and 1+ are simply writ- 
ten in this case using the NRG states of H [see Eqs. (JHJ 
and (|ll|l ] . Below we apply this completeness relation for 
the initial Hamiltonian Ti. 1 , with one slight modification. 
Shifting the to' — to term from Eq. © to Eq. lfTl"|) one 
has the identity 





1 = X~ + 1 + 



(27) 



where 




p 

P 


P P 

1 a m+2 a 
P P 



m+3 



p 
p 


Po 



m — 1 



E 



and 



ei;m| 



^i ,Gi 



FIG. 2: (a) Diagrammatic representation of the recursion re- 
k> e ii m )dis dis(h> e i'i m \ (28) lation of Eq. igOJ for pf^(m). Each box stands for a matrix 

element P,/ ; . [a m +i] (upper row) or its complex conjugate 
(lower row). The state labels U and l[ axe plotted horizon- 
tally, with li to the left and ^ to the right. The state label 
ct m +\ for the m + 1 site is plotted vertically. A connected 
l*"/ line indicates summation over the corresponding index. The 
reduced density matrix p T ^°.(m) has two indices Si and n, 
carried by the external legs on the left. Iterating (a) up to 
m' = N relates the reduced density matrix p T s e . d r°(m) to the 
equilibrium NRG density matrix (A;^; A^|po|fci, A'"). This pro- 
cess is visualized diagrammatically in (b). The sequence of 
connected vertical lines amounts to tracing over the environ- 
ment degrees of freedom {ai}iL m +i- 



T m = E \ k i^u m ){h, 



Here the index ki in Eq. (|29|) runs over all NRG eigen- 
states of iteration m, whether discarded or retained. This 
shift of notation from l m to I m is of purely practical na- 
ture as it allows us to sum freely over all states ki of any 
given iteration to. The symbol X is used to emphasize 
the different complete basis sets of H? and W . 

We are now in position to address the reduced density 
matrix p)j cd (m). Starting from Eq. (|24[1 . we make use of 
the completeness relation of Eq. I|27|l in order to replace 
po with 



(l™+Zjn)po(l-+l?r 



(30) 



is conveniently computed in the course of the NRG it- 
erations as described in Appendix ^ With the matrix 
S(m) at hand, p++(m) is obtained by a simple rotation of 
the "unperturbed" reduced density matrix into the new 
basis: 



Equation (|24|1 is decomposed in this fashion into four 
contributions: 

P™r( m ) = pt,t (irf+Ptr (m)+p7 t ?(Tri)+p-~(m) , (31) 
where pfP (to) equals 

flfrim) = $>, e ; ™K/5oZ£>, e; m) (32) 

e 

(p,p' = ±). Of these four components, only p++(to) can 
be directly related to the "unperturbed" reduced density 
matrix of Eq. I|26|) . To see this, we insert Eq. (|29|l into 
Eq. lO to obtain 

pf,t( m ) =X!E E ( s ' e; m \ k ^ e ^ m > 

x (kl, e' t ; m\p \h, e,; m){ki, e^; m\r, e; m). (33) 

This expression features the overlap matrix elements 
(ki, ef, m\r, e; m), which are diagonal in and independent 
of the environment degrees of freedom: 



p++ 



(to) 



ki,k' 



s i[,s( m )p T kX,kS m ^ Sk ^ r( ^> ■ ( 35 ) 



(ki,ei-,m\r, e;m) = Se^Sk^A 171 ) 



(34) 



The reduced matrix S(m) records the overlap matrix el- 
ements between the NRG eigenstates of TL l m and W m . It 



Note, however, that this transformation is not unitary 
due to basis set reduction at each iteration. Hence, 
pfr( m ) follows directly from the knowledge of p l s c ^ ,0 (m). 

In contrast to pj"+(m), the remaining components 
P+~(to), pj^(m) and pj~(m) cannot be related in 
any simple way to the "unperturbed" density matrix 
// ed '°(TO,). This is readily seen by inserting Eq. (|28|l into 
Eq. (|32|) . which gives rise to overlap matrix elements of 
the form dis( k i, ef, m'\r, e;m) with m! < to. The lat- 
ter matrix elements generally depend on both ei and e, 
which prevents from collapsing the trace over e^ to the 
form that appears in Eq. (|26|l . At present we have no 
feasible way to carry out such a trace. 

Physically, the terms p -1 (to), p '"(to) and p (m) 
encode the way in which high- and low-energy eigenstates 
of TC 1 are coupled within po- For these terms to be im- 
portant, po must contain a significant contribution from 
high-energy states, which means starting from a state 
well removed from thermal equilibrium. In the following 
we restrict attention to the case where po corresponds to 
thermal equilibrium, or more generally to the case where 
Po = InPo^n- 
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2. Starting from thermal equilibrium 



Upon starting from thermal equilibrium, po takes the representation (JSSJ), i.e. (U,ei]m\po 



iteration m. This is a simple consequence of the orthog- 
onality of our basis-set states ji^e^m) and the spectral 

,e';m) = 0, for all 



standard form 



-/97i" 



/Zi, where Zi = Tr{ 



-pn l 



} is the 



unperturbed partition function. In terms of the appro- 
priate NRG basis, it has the spectral representation 4 



Po 



1 



E- 



-0B? 



\k;N)(k;N\ 



Zi = E< 



-0E? 



(36) 



(37) 



These expressions neglect all states discarded at itera- 
tions m < N, exploiting the fact that the discarded states 
have only an exponentially small contribution to po at 
the temperature Tjv^ This NRG approximation usually 
works very well, and can be systematically improved by 
increasing the number of states retained at the conclu- 
sion of each NRG iteration. We denote the latter number 
of states by N s . 

By construction, Eq. (|31)|) obeys the operator identity 
Po = T^Po^n- This guarantees that pf~(m), p~+{m) 
and p~~{m) identically vanish for all m < N. Hence 
p™ d (m) = p++(m) is fully determined by p%f£(m) ac- 
cording to Eq. (|35|) . We now elaborate how p red '°(m) is 
computed recursively from p rod ' (m + 1). Together with 
the "initial" condition 



cw = 5 > 



-/3E* 



(38) 



this provides us with p Ted '°(m) for all m < N . 

Consider an arbitrary m < N . To execute the sum over 
e-i in Eq. I|26|l . we set e, = (a m+ i,e' i ), where e\ encodes 
the N — to + 1 site labels a m +2, • • • , ctjv . In other words, 
e[ is a state variable for the environment R m +x,Ni see 
Fig. [I] Substituting X+ +1 p Xm+i in for p , and using 
the overlap matrix elements 

(si,e i ;m\k' i ,e' i ;m+ 1) = Py.. Si {a m+ i] (39) 

[see Eq. iJBJ], Eq. 1J2T))| is expressed as 

ret 

Xpj^V+l). (40) 

where the sum over ki and fc^ is restricted to the states 
retained at iteration 777 + 1. (For 777 = N — 1, the sum 
runs over all the states of the final NRG iteration). In 
the standard case of a real Hamiltonian W, the matrix 
elements Pn-ijflim+i] are likewise real. Under these cir- 
cumstances one can omit the complex conjugate from 
P£.. r .[a m+ i]. It is also easy to verify that p™ d r°( TO ) van ~ 
ishes if at least one of the states Sj and r% is discarded at 



777 < N . The recursion relation can also be visualized 
diagrammatically depicted in Fig. |3 as brought to our 
attention by Weichselbaum. 

Note that the recursion relation of Eq. (|40|l is not 
confined to thermal equilibrium. It relied solely on the 
fact that po — ^-nPo^-N' which guaranteed that po — 
^■m+iPo^-m+i ^ 0r a ^ rn < N. No additional restriction 
was made to an equilibrium distribution. Nevertheless, 
one is typically interested in the case where the system 
starts from thermal equilibrium. The "unperturbed" re- 
duced density matrix p rcd ' (?77) coincides then with the 
one originally introduced by Hofstetter for the purpose 
of calculating equilibrium spectral functions within the 
NRG^ 



D. TD-NRG algorithm 

After the exposition of the different components of the 
TD-NRG approach, we now turn to the integrated algo- 
rithm. To implement the TD-NRG at temperature T, 
one first selects a chain length N such that T w TV. 
Two simultaneous NRG runs are then performed, one 
for 7i l and another for "Hf . All NRG eigen-energies of 
Uln and Tifn are stored up to the final iteration N. At 
each iteration 777, the overlap matrix S rijr (m) of Eq. I|34|) 
is calculated between all eigenstates 1 7^,777) and \r,m) 
of 7i l m and H^, respectively. Details of the calculation 
are elaborated in Appendix^ This information, as well 
as the product matrices Pi>.i[a m ] for the initial and final 
Hamiltonian, are stored on the hard drive. After the final 
NRG run, the equilibrium density matrix is calculated 
from Eq. H3tj|) using the eigenstates and eigen-energies of 
last NRG iteration for H l N . 

At the conclusion of these steps, the TD-NRG proceeds 
by backward iterations starting from m = N. For each 
backward iteration from 777 to 777 — 1 , the following three 
steps are performed: 



1. The "unperturbed" density matrix p red '°(m — 1) is 
calculated from p red,0 (777) using the product ma- 
trices P;'.,Zi [ot m ] for the initial Hamiltonian TL l m , in 
combination with the recursion relation of Eq. l|4Uf) . 

2. The reduced density matrix p r s cd (m— 1) is computed 
by rotating p red,0 (777 — 1) to the basis of the final 
Hamiltonian using the overlap matrix S(m— 1) and 
Eq. (HJ). 



3. Using Eq. I|25fl . the contribution of iteration 777 to 
0(tj) is evaluated simultaneously for all times tj of 
interest. Here the only matrix elements of p^^m) 
and 0™ s to contribute are those where at least one 
of the states s and r is discarded at iteration 777. 
At the conclusion of this step, p Tcd (m) is deleted to 
reduce the memory load. 



These steps are repeated until m = m m ; n is reached, 
below which no state has been discarded. 

It is easy to see that if Ji* — Ti l , i.e., the Hamiltonian 
is left unchanged, then 0(t) so obtained coincides with 
the thermodynamic average of O for all t. Indeed, the 
overlap matrix S(m) is diagonal in this case, which means 
that p red (m) and p red '°(m) are the same. Consequently, 
p l s C r(m) has nonzero matrix elements only if the states r 
and s are both retained, leaving only the iteration m — 
N to contribute to Eq. (|2"5)) . Since p Icd (N) coincides 
with the equilibrium density matrix which is diagonal 
in energy, one recovers the thermodynamic average of O 
independent of t. 

A word is in order at this point about the nonequi- 
librium NRG approach of Costipi& and its relation to 
the TD-NRG. Costi had focused on the calculation of 
nonequilibrium spectral functions. In contrast to the TD- 
NRG, he settled with a single NRG shell m to evaluate 
the spectral function at frequency to ~ D m . In practice 
this meant replacing the reduced density matrix p™^ (m) 
with the full equilibrium density matrix at iteration m. 
After rotating the latter according to Eq. 1(35(1 . only a 
single Wilson shell contributes to a given uu of the spec- 
tral representation of 0(t). Being well aware that a full 
multiple-shell evaluation is ultimately required for the 
correct description of nonequilibrium dynamics^ Costi 
carefully applied the single-shell approach only to inter- 
mediate frequencies. This restricted his ability to track 
the real-time dependence of 0(t). 

In the TD-NRG we overcame these difficulties by iden- 
tifying an appropriate basis set for the Hilbert space of 
the iV-site chain, and by the resummation procedure that 
has led to Eq. ((25(1 . The significance of these steps is best 
reflected in the fact that 0(t — ► 0+) of Eq. (|25j) exactly 
coincides for any local operator O with its NRG ther- 
modynamic average with respect to 7i l N - While physi- 
cally clear, this statement is highly nontrivial from the 
standpoint of the TD-NRG, as all energy scales con- 
tribute to the summation of Eq. ((25(1 . Moreover, the 
limit 0(t — > oo) stems from the degenerate terms with 
E" 1 — E" 1 from all iterations m in Eq. ((25(1 . In frequency 
domain, these terms give rise to an addition 5 (us) contri- 
bution that comes out naturally in the TD-NRG, but is 
absent in Costi's approach. The latter is of crucial impor- 
tance for obtaining the correct asymptotic long-time be- 
havior of operators with non-vanishing expectation val- 
ues. 



E. Towards restoring the continuum limit 

So far, we focused on accurately calculating the time 
evolution of any local observable O on a finite-length 
chain. While this might be a reasonable representation of 
mesoscopic systems characterized by a finite level spac- 
ing, our approach is geared toward the description of a 
continuous bath. Relaxation in such a macroscopically 
large system cannot be fully described by a finite Wilson 



chain. Indeed, it was already pointed out in the context 
of the equilibrium NRG™ that certain thermodynami- 
cal properties undergo unphysical oscillations as a result 
of the discretization of the continuous bath. To circum- 
vent this problem, Oliveira and coworkers^ introduced 
a ^-dependent logarithmic discretization of the continu- 
ous bath according to [1, A~ z , A~ z ~\ • • • , A^ z ^ n ~ 1 ,- ■ ■]. 
Wilson's original discretization corresponds in this notion 
to z = 1. As shown by Oliveira et a/., the unphysical 
oscillations can be removed by integrating expectation 
values with respect to < z < 1. This has the effect 
of mimicking a continuous bath using a single adjustable 
parameter. 

We apply the same technique to improve on the compu- 
tation of time-dependent quantities. To partially mimic 
the relaxation in an infinite-size system, we calculate the 
time evolution of Eq. (|25|l for each value of zt — i/N z , 
i = 1, • • • , N z , and average over the different z^s. Here 
N z is the total number of z-values considered. Since the 
NRG oscillations are primarily associated with sin(27rz) 
and cos(27rz) terms on the interval < z < lr& N z 
should be chosen in multiples of 4. This leads to optimal 
cancellation of oscillations. As shown below (see Fig.[3J), 
usage of the z-trick can greatly improve on the long-time 
behavior of the TD-NRG. 

Even with the z-trick at hand, the evaluation of Oit) 
boils down to summation over a finite number of oscilla- 
tory terms of the form e l ( Er ~ E ° )'. One way to simulate 
a continuous spectrum is to broaden the NRG levels, as is 
done in the calculation of equilibrium spectral functions. 
This requires, however, extra care in the nonequilibrium 
case. As mentioned above, the limit 0(t — » oo) origi- 
nates from the time-independent terms with E™ = E'" 1 
in Eq. 1)25(1 . These terms must remain in tact in order 
to correctly describe the long-time behavior. Therefore, 
one must separate the time-independent terms from the 
oscillatory ones. Damping should only enter the latter 
terms in Eq. J2U). 

Focusing on the oscillatory terms, we replace each en- 
ergy E™ with a Lorentzian broadening according to 



„±i_E, m t 



dE 



r, 



„±iEt 



,±i_B, m t-r m t 



tt (E-E^f+Tl 



(41) 

Here T m = a c iD m with ad a constant of order one is a 
scale-dependent broadening, reflecting the characteristic 
energy resolution of the mth NRG iteration. In this fash- 
ion, each of the terms E™ ^ E™ in Eq. (|23|) is modified 
according to 



i(E™-E™)t 



a^iK 



-E™)t-n d D m t 



(42) 



while the terms with E™ = E™ are left in tact. Below we 
present results with and without the additional damping 
factor ad- 
It should be noted that the Lorentzian broadening of 
Eq. I|41|) differs from the customary log- normal form used 
for equilibrium spectral functions. This choice is moti- 
vated by physical considerations, as it produces a simple 
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exponential decay in time. We also stress that ay comes 
to simulate the continuous spectrum of the bath. Should 
Ti possess eigenstates which are pure eigenstates of the 
impurity part of the Hamiltonian 7ij mp and do not couple 
to the bath (as is the case in certain parameter regimes 
of electron-transfer models^), ay must be set to zero. 
Else, unphysical damping is introduced. 



F. Discussion of the TD-NRG approach 




analytic 


al result 


A=1.3 




"-■ A=1.3 2 




— A=1.3 3 




— A=1.3 4 




A=1.3 4 


N z =l 



t*r 



For clarity, we summarize the key ingredients as well 
as the assumptions made in the TD-NRG approach in 
this section. We have to divide the different aspects of 
conceptual and technical nature. 

As mentioned in the previous section, one of the most 
important conceptual points to bear in mind is the mim- 
icking of a continuum of bath degrees of freedom by a 
selection of discrete states. Therefore, we expect al- 
ways some oscillatory contributions in the time evolu- 
tion of 0(t): the shorter the chain length, the smaller 
the amount of bath degrees of freedom the stronger the 
oscillations. This property is inherent to any discrete rep- 
resentation of a continuum and not a shortcoming specific 
to our method. In fact, the mapping error is controlled 
by the NRG discretization parameter A and the chain 
length becomes infinite for a given temperature T for 
A — ► 1+. In addition, we could show that these oscilla- 
tions are strongly suppressed by the z-trick described in 
the previous section also mimicking a bath continuum. 

Another conceptual point concerns our expectation of 
the time evolution. If we would include all possible physi- 
cal interactions, have all baths for energy and particle ex- 
change coupled to the system, wait an infinitely amount 
of time, we expect that the system equilibrates to the 
new thermodynamic equilibrium governed by Ti* inde- 
pendent of the initial conditions. This is required by the 
basic assumption of equilibrium thermodynamics. How- 
ever, this is in general not the case. Imagine an isolated 
spin in metallic host in strong magnetic field, as we will 
investigate later in the section on the Kondo model. If we 
artificially decouple this spin from the environment and 
switch of the external magnetic field, the local spin will 
remain its polarized steady state rather than relaxing to 
the new thermodynamics state since the local spin re- 
mains a conserved quantity. The local expectation value 
can only evolve towards the new thermodynamics equilib- 
rium, if H? provides an energy, magnetization or particle 
exchange mechanism. 

Imagine, A7i is controlled by an external time- 
dependent field f{t) which couples to a conserved quan- 
tity Q of the total system, i. e. [Q,W] — and AH = 
f(t)Q. Then all eigenstates of TC l remain eigenstates to 
H.f for arbitrary times. Only the eigen-energies obtain 
a time-dependent shift £ m (t) = E l m + f(t)q, where q is 
the eigenvalue of Q of state \m). It is obvious that there 
will be no time dependence of any operator whose matrix 
elements are energy-independent such as spin or charge. 



FIG. 3: Comparison of the time-dependent occupancy nd{t) 
of the resonant-level model, obtain analytically by Eq. 1521 
and numerically by the TD-NRG through the evaluation of 
Eq. $£%. A sudden change of E d from E° d = to E d = -2r 
is considered without an accompanying change of the hy- 
bridization strength: To = V\ = I\ Different values of 
A = 1.3, 1.3 2 , 1.3 3 , 1.3 4 are used. The number of NRG iter- 
ations (i.e., the Wilson chain lengths) have been adjusted so 
that all curves are calculated at approximately the same tem- 
perature Tn- Explicitly, in going from A = 1.3 to A = 1.3 4 , 
Tjv/r equals 0.00171,0.00176,0.00183,0.00193, correspond- 
ing toN = 96, 48, 36, 24 NRG iterations. The remaining NRG 
parameters are as follows: N 3 - 1000, N z = 16, D/T = 500, 
and a d = (i.e., no additional damping). 



Another non-trivial difference to the assumption of 
thermodynamics arises from the fact that we describe 
the time-evolution of a closed system. This is in contrast 
to the Keldysh approach 3 which introduces - usually in 
a very subtle - an infinitely small relaxation rate much 
smaller than any energy scale in the problem to ensure 
the asymptotic approach of the thermodynamical state. 
Sometimes, even explicitly the existence of an additional 
thermodynamic bath is made. In our description, how- 
ever, the total energy of the system remains a constant 
and is given by 

(H)(t>0) = Tr{pit)H(t)}=Tr{p n f } (43) 

since the time evolution operator exp(— ihi^t) commutes 
with Ji^ . In contrary, the total energy of the thermody- 



namical equilibrium is given by 



in 3 ) = Tt{p f {p)H J ) 



f^ 



(44) 



where £/(/?) = exp(-/3H / )/Tr{exp(-/37i: / )}. A sudden 
change in the Hamiltonian of closed finite size system 
will always result in an effective heating, where the ef- 
fective temperature /3 e // can be obtained by solving the 
implicitly equation 



Tr{p (P)H f } = Tr{p f (0 eff )H'} 



(45) 



If the energy of the bath states of H l 'f, however, 
are continuously distributed and their energies are not 
changed by ATi, the temperature will not be changed in 
the thermodynamic limit. 
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III. FERMIONIC BENCHMARK: RESONANT 
LEVEL MODEL 

A. Definition of the model 



equal-time lesser Green function G d (t,t') = (eft (t')d(t)) . 
Using the Keldysh technique (i*2i we derived an exact an- 
alytical expression for rid(t) in the wide-band limit, for a 
step-like change of parameters at t = 0: 



The resonant-level model is perhaps the simplest exam- 
ple for a quantum-impurity system. A single Fermionic 
level d) with energy Ed is coupled by hybridization to a 
bath of spinless Fermions: 



H = ^e k c\c k + E d {t)d)d 

k 



(46) 



Here c k creates an electron with momentum k in an s- 
wave state about the position of the level (taken to be the 
origin). To make connection to the TD-NRG approach, 
we consider a step-like change in Ed(t) and V(t): Ed(t) = 
E° d 0{-t) + E d 9(t) and V{t) = V 8(-t) + V x 6{t). Since 
the model is bilinear in Fermionic operators, equilibrium 
properties (i.e., for E d = E d = Ed and Vo = V\ = V) can 
be calculated exactly using, for instance, the local Green 
function Gd(z): 



G(z) 



z-E d - A(z) ' 



A w-E 



JVf_ 
z- £fe 



(47) 



(48) 



In the wide-band limit, T = Im{A(— i0 + )} defines the 
width of the level due to the coupling to the Fermionic 
bath. 



B. Analytical Keldysh solution for nonequilibrium 

Out of equilibrium, the occupancy of the level, n<j(£) = 
(d'(t)d(t)), follows directly from the knowledge of the 



n d {t)= PF f(e)\A(e,t)\Ue 



(49) 



with 



A(e,t) 



dr 



y/ T \ e -i*r-if*d(;lE d (l;)+ir(Z)] 



Vi 



Ti + i(E d - e) 

Vo 



-(ie+ri)t 



Vi 



r + i(E° d - e) T x + i(E d - e) 



(50) 



Here, p F — X)fc^( £ fe) ^ s tne density of states of the 
Fermionic bath, /(e) is the Fermi-Dirac distribution func- 
tion, r(£) equals Trp F \V(t)\ 2 , and Vi = irp F \Vi\ 2 (i — 
0,1). 

Besides the wide-band limit, Eqs. I|49l) and l|50|l require 
that | Vb| > 0, which comes to ensure that the initially de- 
coupled level for t — > — oo has decayed to its equilibrium 
state for — oo < t < 0. Since any infinitesimal value of Vo 
fulfills this requirement, the limit V"o — > can be viewed 
as contained in Eqs. (|49() and i|50|) . 

For T — v 0, Eqs. (|4*9"|l and JSUJ) can be evaluated in 
closed analytic form. Introducing the auxiliary function 



F(x,y;t) = -e tv+itx Ei(ty + itx) 

+2iri6(-x)6(-ty)sign(t)e ty+ltx , (51) 



where E\(z) is the exponential integral function, we ob- 
tain 



n d (t >0,T = 0) = (l + e- 2rit )n d + e- 2rit n d 

-2T t t 



-2syToTi Re 

7T 

-2 Sl /r ri Re 



E d-E d -i(r + r 1/ 



e ~iE d t 



-rit 



[ln^S-iToJ-lnC-Ei+iTi)] 



+2- 



-Im 



I 



[F(E° d ,T ;t)~F(E d ,-r i; t)] 



(52) 



Here n d = \ — — arctan(_E^/Fi) with i = 0, 1 is the equilibrium occupancy of the level for the corresponding 
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model parameters, and s — sign(ViVo)- 

We note that the exact result does not exhibit a sim- 
ple exponential decay to the new equilibrium occupancy 
n d . It is actually governed by two relaxation rates, Fi 
and 2T\. Moreover, the last and second- to-last terms 
in Eq. (|52[1 contain an oscillatory factor exp(— iE^t), 
responsible for Rabbi-type oscillations that visible for 
|2?i|/ri > 1 at short time scales. In a previous publi- 
cation^ we have used this rather elaborate formula and 
its extension to finite temperatures to benchmark the 
TD-NRG. Here we complete the discussion by focusing 
on the role of the NRG parameters A and N z . 

Figure presents a comparison between the exact ana- 
lytic result of Eq. I|52l) and the time-dependent occupancy 
n<i(t) obtained by the TD-NRG for different values of A 
and N z . The temperature is roughly kept fixed in all 
curves by adjusting the length of the NRG chain. While 
the curve corresponding to the smallest value of A and, 
thus, to the longest chain used (N = 96) shows excellent 
agreement with the exact analytical result on all time 
scales, the curves for the larger values of A show good 
agreement only at shorter times, t T < 2. Deviations 
from the exact result develop at longer times, which are 
characterized by oscillations about a value of rid that is 
reduced as compared to the new thermodynamic average. 
We emphasize that none of the TD-NRG curves in Fig. El 
involved an extra damping factor ay > 0. 

The above behavior is not surprising since, as discussed 
in Sec. Ill Fl the shorter the NRG chain, the worse it rep- 
resents a continuous bath. A different way of viewing the 
matter is through the conservation of energy and parti- 
cles in the system for t > 0. In a macroscopic bath, 
the system relaxes to the new equilibrium state by redis- 
tributing the excess particles among the different lattice 
sites. Since the excess number of particles is of order 
one, each lattice site (the level included) acquires an in- 
finitesimal shift to its occupancy, which vanishes in the 
thermodynamic limit. For a finite system, the shift in 
occupancy remains finite. We therefore conclude that 
the oscillations depicted in Fig. [3] are real in a finite-size 
system. However, they are unphysical for a continuous 
bath, as demonstrated by the exact analytical solution. 
We pinpoint this "inaccuracy" to be of conceptual na- 
ture, controllable by the limit A — ► 1 + . 

The magnitude of the finite-size oscillations is greatly 
reduced, though, by resorting to Oliveira's z-trick. To 
illustrate this point, we have added one extra curve for 
A = 1.3 4 = 2.8561 without any z-averaging. Enhanced 
by the short chain length, the finite-size features are no- 
tably more noisy and have a larger amplitude. The effect 
of averaging over z is to greatly smoothing and reduce 
these finite-size structures. Hence averaging over z is vi- 
tal for describing the long-time behavior, particularly if 
no extra damping is introduced. 

The short time behavior, however, is very accurately 
reproduced defying prejudice against the NRG which al- 
legedly is insufficient to describe the high energy physics 
correctly needed for the short time behavior. Even if 



the statement would be true, the NRG represents the 
high energy physics incorrectly for calculating thermody- 
namics properties, which is not - see Wilson^ - it would 
have no influence on the time evolution calculated with 
the TD-NRG for a very simple reason. The high energy 
states are generated at the first iterations. For energy 
scales much larger than the relevant scales of the quan- 
tum impurity, the NRG flow of TV 1 and Jt* are more or 
less identical which leads to a almost diagonal overlap 
matrix S r _ s (m) for these iterations. Since the reduced 
density matrix p™ ' k ) (m) is only non-zero for retained 
states k, k' there will be no contribution to the time evo- 
lution when evaluating (|25|l . Only when the flow be- 
tween the two Hamiltonians starts to deviate, contribu- 
tions to the expectation value 0(t) are generated. This 
may be seen clearly in the inset of figure 2a of our previ- 
ous paper— which shows the iteration dependent evolu- 
tion of rid(t — > + ). 



IV. BOSONIC BENCHMARK: DECOHERENCE 
IN THE SPIN-BOSON MODEL 



The spin-boson model2£*22 

A 



H 



&X + -er 2 + ^^afa» 



tEm«! 



(53) 



is perceived as the simplest model for the understanding 
of dissipation and decoherence in quantum systems. A 
two-level system, represented by a spin, is coupled to the 
displacement operator of a continuum of bosonic degrees 
of freedom with the coupling constant Aj to each bosonic 
mode. This continuum causes decoherence, dissipation 
and possible spin decay. The effect of this dissipative 
environment is fully captured by the coupling function 
J(w) 



J(u>) — Tr} y A?fl(o> — ujj) . 



(54) 



Since the low-energy part of the spectrum dominates the 
response of the spin, the high-energy details of J(u>) can 
be neglected: this function is usually approximated by a 
power-law behavior—: J[ui) oc ui s , < u < w c . In the 
literature, one distinguishes between the "ohmic" case, 
,s = l, corresponding to a classical model with a dissipa- 
tive term proportional to the velocity, the "sub-ohmic" 
case < s < 1 and the "super-ohmic" case s > 1. In the 
ohmic regime, the spin-boson model can be mapped onto 
the anisotropic Kondo model in the wide-band limit^ 
with A/uj c — pfJ± and a = (1 — pfJz) 2 - 

This model plays also an important role for under- 
standing charge transfer reactions in chemistry^ where 
the spin states would represent the initial and final state 
of the transfer reaction and the bosonic continuum the 
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FIG. 4: Decoherence of the off-diagonal density matrix el- 
ement poi of the local reduced density matrix in the spin- 
boson model as function of time, by measuring the expec- 
tation value s x (t) = poi- The solid lines show the TD- 
NRG results, the dashed line with the same color (online) 
the exact analytical result given by 1551 . NRG parameters: 
N s — 150, Nu e r = 14, Nb — 8 (number of bosons added a each 



iteration 



N z = 
0,a d 



16, A 
= 0.1. 



V2.T = 0.0078, a = 0.1, w c = 



density fluctuations of solvents or molecular vibrations. 
It has been argue d 30 ' 31 that the model serves a quantum 
mechanical foundation of the classical Marcus theory^ 
for electron transfer reaction. 

In addition, the model has been discussed in the con- 
text of decoherence of qubits in quantum computing 3334 . 
Setting A = 0, a z is a conserved quantity, and, therefore, 
the bosonic environment can be traced out analytically. 
If we prepare an initial density matrix containing locally 
only the pure state \s) = (|0) + |1))/a/2, it was show: 
that off-diagonal part of the reduced density matrix 



34 



pio(t) = (l|Tr BoS o„{p(i)}|0) 



(55) 



can be written as pw(t) = e r ^pio(0). Here, r(i) is 
given by the exact expression 



T(t) = 



du: J{u>) coth 



to \ 1 

27V ~~ 



cos(urf) 



(56) 



for the temperature t£&2^ 

Recently, Bulla and collaborators introduced an ex- 
tension of the NRG to bosonic baths^ 1 ?' 3 ^ These authors 
investigated the different quantum critical points of the 
model as a function of the bath exponent s and the 
coupling strength a. We combined their bosonic NRG 
method and our TD-NRG approach to show explicitly 
that our method is very accurate for short as well as long 
times up to t * T s=s O(l) for bosonic baths as well. 

In order to start with the same boundary conditions 
as the exact result was derived for, we set A = 100 in 
7i l . This guarantees that equilibrium density operator 
consists of the disentangled operator product of \s)(s\ 
and the equilibrium density matrix of the bosonic bath 
for t < 0. For t > 0, A is set to zero and an entangled 
density matrix is evolving with time. The off-diagonal 



matrix element pio corresponds to the expectation value 
of s x (t) obtained by the bosonic TD-NRG. The coupling 
function J(ui) 



J{w) — 2-iraLo]. s uj s for < u) < uj c 



(57) 



enters the analytical as well as the numerical calcula- 
tion. We have obtained Poi(*)/Poi(0) for fixed a = 0.1 
and different exponents s. The comparison between the 
TD-NRG and the exact analytical result given by Eq. 
(|55|1 is depicted in Fig. 0] for the super-ohmic, the ohmic 
and the sub-ohmic case. We observe excellent agree- 
ment between the exact analytical result and the pre- 
dictions of the TD-NRG p 4 ^ The small oscillations in the 
quantum regime observed in all curves in the time range 
0.01 < t *T < 1 are no artifacts of the method but real 
features. By comparison with the exact analytical result, 
they have been traced to the usage of the hard frequency 
cutoff in J(lu) at lu — lu c . Using a soft cutoff function 
J(lu) cx lu s cxp(— lo/uj c ) would produces smoother and 
slightly shifted curves. 34 

We noted that due to the larger numbers of degrees of 
freedom added though each chain link, the time evolu- 
tion for a bosonic bath shows a higher accuracy and less 
dependence on A than for fermionic baths as long as one 
does not leave the range of validity of the bosonic chain 
NRG^ 



V. THE KONDO MODEL: ANALYTICAL 
CONSIDERATIONS 

A. Definition of the model 

The Kondo model comprises of a local spin Si mp in- 
teracting with a spin degenerate conduction band via a 
local anisotropic Heisenberg term 

Hr = Jz 2_j ac R=0.<j C R=0,<? S imp (58) 

a 

+ J ± \ C R=0,1 C R=°-l S vmp + C R=0,l C R=W S tmp) 



where 



J 

-R=0,a 



creates a conduction electron at the posi- 



tion of the spin R = 0. For J z = J±, we obtain the 
usual SU(2) symmetric interaction. We also allow for a 
time-dependency of the coupling constants J z and Jj_. 
The energy of the local levels are splitted in an external 
magnetic field H(t) by the Zeeman energy 



H,, 



-lsH{t)S v , 



(59) 



whose transversal parts can also be interpreted as hop- 
ping terms of the two-level system modelled by the lo- 
cal spin 1/2. The local conduction electron creation op- 
erators are expanded in s-wave eigenmodes of the non- 
interacting conduction electron band 



"R=0,a 



(60) 
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where p(e) denotes the conduction electron density of 
states^ Note that the anti-commutator of the local con- 
duction electron operator is given by {c* R=0 a , cr=o, CT '} = 
8 aa i . The conduction electron Hamiltonian is given by 



H, 



E 



de e c\ a c €a 



(61) 



where c\ a (c e( ,) creates (annihilates) electron in an s-wave 
state with spin a, energy e. All other angular momen- 
tum do not interact with the local spin^ The integration 
runs from the lower to the upper band edge, Di ow to D up . 
Throughout the paper, we assume a symmetric band with 



a relativistic dispersion, i. e. D 



Ion: 



-D 



up 



-D and 

Pf = p(e) = const = 1/(2.0). Then, the total Hamilto- 
nian under consideration is given by 



H(t) = H c + H K (t) + H loc {t) 



(62) 



Why do we neglect the spin polarization of the con- 
duction electrons? As we discussed in section 111 Fl an 
external field coupling to total z-component of the spin 
yields a time independent local spin expectation. Since 
the external magnetic field, however, couples to the total 
magnetization 



Mtot = lsS a + j b St ot 

= (js-7b)s s + lb S tot 



(63) 



where s z is the local spin operator and Sl ot the total 
spin of the conduction electrons, the relevant Zeeman 
energy for the real time dynamics is given by H zee = 
— ("f s — 7(,) s s H. In this case, the local spin can relax if 
the gyromagnetic ratios of local spin 7 S and conduction 
electrons 7b differ. Then, j s would have to be replace by 
7s = ils-Jb) inEq. {52). 

In real materials, however, the spin-orbit scattering of 
the conduction electrons on distributed impurities in the 
metallic host furbishes another relaxation mechanism, as 
pointed out by Langreth and Wilkins 1 more than 30 years 
ago. We might argue in the following way. If this relax- 
ation processes for the conduction electron magnetization 
is much faster than the relaxation process of the impurity 
spin due to coupling J to the conduction electrons, we 
can consider the conduction electron band as demagne- 
tized on the time scale of the impurity relaxation process. 
In this case, neglecting the spin-polarization of the con- 
duction electrons are justified. If for some reason, the 
spin-lattice relaxation processes is very slow compared 
to the impurity spin, the impurity spin-relaxation pro- 
cess is dominated by the spin-lattice relaxation time in 
which case we have to substituted Y s = (7 S — 7&) or rescale 
the absolute values of the physical magnetic field. 

It turns out, that the long-time behavior of impurity 
spin relaxation is dominated by the Kondo time tx in 
the Kondo regime, given by the reciprocal characteristic 
thermodynamic energy scale tic = 1/Tk, and by the re- 
ciprocal temperature 1/T for the ferromagnetic regime 



of the model. As long as the spin-lattice relaxation time 
is finite, for sufficiently low temperature T and Kondo 
temperature Tr , neglect of the conduction electron spin- 
polarization is justified. However, the short time behav- 
ior of the spin relaxation has admittedly purely academ- 
ical value if t sp > 1/J. 

We will distinguish between two different starting po- 
sitions. Either we leave the coupling constant unaltered 
and just switch of the magnetic field at t = 0, which is 
the physically more relevant condition. Alternatively, we 
start from an initially decoupled system, the local mo- 
ment fixed point of the Hamiltonian, and switch on the 
Kondo couplings J z and J± as well as changing the ap- 
plied magnetic field at time t — 0. It was previously 
pointed outi^& that this setup is closer to the thermo- 
dynamics since the time t might play a similar role as 
/?. 



B. Analytical results 

Wilson invented the numerical renormalization group 4 
to solve the thermodynamics of this model very accu- 
rately for all parameter regimes. The model has two fixed 
points for the antiferromagnetic regime, the unstable lo- 
cal moment fixed point (Jz = 0, J±_ = 0) corresponding 
to P = 1/T — ► and the stable strong coupling fixed 
point ( J z — > oo, Jj_ — > oo) for (3 — * oo caused by the 
infrared divergence in the perturbation theory in the ab- 
sence of an external magnetic field. For the ferromagnetic 
regime, i.e. J z < and \J Z \ < \J±, the stable fixed point 
is given by (J z (oo) < oo, J± = 0). Due to the ferromag- 
netic polarization of the conduction band in the vicin- 
ity of the impurity, the spin-flip scattering is more and 
more suppressed for decreasing temperature: the model 
remains perturbative in this regime. 

The analytical "poor man's scaling" second-order 
renormalization group treatment^ prior to the NRG pre- 
dicts 



J z — J± = const. 



(64) 



The renormalization of isotropic coupling J(D') as of the 
running cutoff D' 



J(D') 



J(D) 



1 - J\n{D/D') 



(65) 



provides some analytical inside for the fixed point 
structure even though the treatment breaks down for 
Jla(D/D') -»■ 1. 

For the isotropic ferromagnetic Kondo model, however, 
J{T) remains finite and vanishes for T — > 0. Since the 
model approaches the local moment fixed point in this 
case, a spin-polarized local state remains spin-polarized 
after switching of the magnetic field at T — without 
any additional relaxation mechanism present. At finite 
temperature, however, a thermal energy of T provided 
fluctuations in the bath in addition to the very small but 
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finite transversal spin flip term J±(T) = J(T). There- 
fore, we expect that the FM Kondo model has a charac- 
teristic time scale t F M = \2p F J{T)T\~ x which is actually 
found in the calculations (fi = 1). Through the paper we 
measure the temperature in energy units (fcg = 1.) 



C. Short time behavior and relaxation towards the 
Kondo strong coupling fixed point 

If we consider a decoupled impurity for times t < 0, we 
can calculate analytically the short time response using 
a perturbation expansion. For J z = J z (t < 0) = and 
Jj_ = J±(t < 0) = 0, the equilibrium density operator is 
given by exp(~f3HQ)/Z where H® = H c +Hi oc . Switching 
on the coupling J\ = J z (t > 0),J| = J±(t > 0) at t = 



leads to the time evolution of spin operator 



Sz(t) 



iHl 



S* 



-iHt 



(66) 



where H = H c + Hk , if the magnetic field is switched off 
for t > 0. In the interaction picture, the spin operator 



Sl{t) = e- lH ^S z {t)e 



obeys the equation of motion 



m c t 



dSJjt) 
dt 

which is integrated to 



= i[H K (t),Si(t)} 



(67) 



(68) 



Jo Jo Jo 

; / dn f 1 dr 2 fVsIi^TiMi^Mi^Ts),^)]]] 

Jo Jo Jo 



(69) 



+z° 



where -ff^(t) = e~ tHct HKe lHct . Neglecting the last term 
yields a perturbative result correct up to order 0(J 2 ). 
Since S z commutes with the S z term in Hk, its short 
time behavior is solemnly governed by J±. The renor- 
malization of J± by J z sets in at third order, well known 
for the Kondo problem^ Alternatively, the time evo- 
lution of (S z (t)) could have also been calculated using 
Keldysh techniques^ All expectation values have to be 
evaluated with equilibrium density operator po- Since 
we are only interest in an analytical results for short 
times, the Keldysh approach does not have any advan- 
tages: it also fails to capture the Kondo physics at long 
time scales i 

There will be no contribution linear in 0(J Z , JjJ since 
(S + ) = (S~) = 0. Evaluation of the two commutators 

needed for the second order contribution ASi in l|69() 
yields 

/>oo />OC 

A^ 2 ) = (S z )(2tJi) 2 / de / de'p(e)p(e') 



— oo J — oo 



x/(e')(l-/(e)) 

xSRe 



it 



(e-e') 



A2 



(70) 



For T — > and a constant density of states pf = 
1/(2D), we can solve the integrals analytically and obtain 



(S z (t)) - (S z 



ASi 2 



(71) 



r 



where we have defined the function G(x) by its series 
expansion 



G(x) = £ 



( _ 1)J+ i 



/=i 



(21)121(21-1) 



(72) 



(S z )(2 PF J\_) 2 [G(2Dt)-2G(Dt)] 



The result is exact up to 0(J^_): the short time evolution 
is independent of J z , consistent with the physical picture 
that a change of the spin state can only by mediated 
by the spin-flip term of Kondo interaction. J z enters at 
first in the third order of the perturbation expansion, 
corrections which are neglected here. 



VI. THE KONDO MODEL: NUMERICAL 
RESULTS OF THE TD-NRG 

A. Antiferromagnetic regime 

1. Short and long time behavior for T — » 

In order to make connection between the TD-NRG 
and the analytical calculation for the short time behavior 
presented in previous section, we start with a decoupled 
impurity in a finite magnetic field along the z-axis, i.e. 
J° = J\ = 0,H z /D = 0.1. At t = 0, we switch on the 
anisotropic Kondo interaction and simultaneously switch 
off the magnetic field. At infinitely long times, we expect 
that the system will relax into the strong coupling Kondo 
fixed point for T — > and t — > oo, independently of the 
initial conditions. 



16 






J 


=0.15 - 


-. J-= 


-0.1 


0.44 


--■ I, 
' .... J z 
.... i y 


=0.1 - 
=0.05 . 

=0.0 

... 1 




=0.15 (analytic) 
=0>15(ana) 0(t 2 ) 

i .'" ' 




t*D 



0.5 

0.4 

"^0.3 

0.2 

0.1 


10 



: 


— 2pJ z =0.15 
-- 2pJ z =0.1 

— - 2pJ z =0.05 
.-. 2pJ ? =0.0 
-- 2pJ z =-0.1 

— Flow Equation 


— '^£i 


- 


>^_,,J 



t*T„ 



FIG. 5: Spin-expectation value for the anisotropic Kondo 
model in the antiferromagnetic regime vs temperature for dif- 
ferent values of J z /D = -0.1,0.0.05,0.1,0.12,0.15 and fixed 
J±/D = 0.15. The upper panel (a) shows the short time 
behavior and the analytic fit using Eq. I|71^ (magenta curve) 
and its the second order contribution in t. The panel (b) 
displays the long time behavior of the same data on a loga- 
rithmic time scale while in (c) the data of (b) isplotted vs 
t/tx ■ We added the universal curve from Ref. [3JJ where we 
rescaled their definition of Tk by the factor 1/1.67 to match 
ours. NRG parameters: A = 1.5; N s = 1000, a = 0.1, N z = 
16,J° = Jl = 0, H°/D = 0.1, Na*. = 130,T = 3.7*iO- 12 . 



In Fig. |SJ the time evolution of S z (£) is shown for a 
fixed spin flip rate J±/D = 0.15 and different values of 
J z /D. S z (t) is independent of J z as expected for times 
t*D < 1, clearly visible in the upper panel (a). Only spin- 
flip processes can alter the spin polarization. At short 
times, the Kondo correlations are absent. Only the bare 
coupling constants enter the analytical result calculated 
for the comparison with the presented TD-NRG curves. 
The magenta (online) curve in (a) displays the second 
contribution in J± given by Eq. I|71|l which contains ar- 
bitrary high orders in t through the series expansion of 
G{x). Since the analytical formula is independent of J z , 
it must correspond to the TD-NRG curve for J z = 0. 
We find an excellent agreement between the analytical 
and the TD-NRG result confirming the accuracy of our 
method at short and intermediate time scales. These 
findings show that the ultra-short time response is purely 
perturbative, since the strong correlations develop only at 
low energies and, therefore, influence only the long time 
behavior. The J z > curves lie below, the curve with a 
ferromagnetic J z /D — —0.1 above the analytical results. 



Jz/D 


T K 


t K = 1/T K 


0.15 


0.000537 


1862.3 


0.10 


0.000203941 


4903.4 


0.05 


5.49093 * 10 -5 


18211 


0.0 


8.00952 * 10~ 6 


124850 


-0.1 


2.31866 * 10~ 10 


4.3128 * 10 9 



TABLE I: Kondo temperature as function of J z for fixed J± — 
0.15 used for rescaling the time axis in Fig. [^Ji). Note that 
the Kondo temperature Tk, obtained by the thermodynamic 
condition^ A(S z )(T K ) = 0.07, has a A dependence. NRG 
parameter A = 1.5, N s - 1000. 



Starting with order 0(J 3 ), J z terms contribute to S z (t) 
by renormalizing J± . The positive J z leads to an increase 
of Jj_, while a ferromagnetic J z yields a reduction of J± 
consistent with our findings in Fig. (Sp) of an increasing 
relaxation time scale with decreasing J z . The long-time 
behavior of S z (t) is displayed in Fig. [SJj) ■ Note that the 
TD-NRG accurately predicts the relevant time scales of 
the spin relaxation for arbitrary long time scales, here 
over 10 orders of magnitude in units of the reciprocal 
band width 1/D as long as t * D < 1. 

In order to identify the long time scale, we rescale the 
data with the thermodynamic Kondo temperature Tk 
which is obtained from the temperature dependency of 
the effective local moment A{S 2 Z ){T K ) = 0.07. Excellent 
scaling is found for the long time tails. Instead of using 
the thermodynamic Kondo temperature, we could have 
defined an effective time scale ti ow by which (S z }(ti ow ) — 
0.1 is reduced to 20% of the starting value S z (0) = 0.5. 
It turns out that ti ow * Tk — const w 0.76 independent 
for J z within the numerical accuracy. Therefore, the two 
time scales are equivalent, and the Kondo time tK = 
1/Tk governs the long time behavior of S z (t). 

From this discussion, it is apparent that there will be 
no universal master curve for S z (t) valid for all times 
and solemnly parametrized by a universal time scale ti ow . 
Only at very long times t ^> tiow, the curves tent to 
approach a universal curve S{t/ti ow ). The panel (c) of 
Fig. presents how the value of J z determines this ap- 
proach to this master curve. We also added the universal 
curve for the weak coupling isotropic Kondo model, i.e. 
PfJ <C 1. Lobaskin and Kehrein have calculated this 
curve 59 using Wegner's flow equation 40 and expanding 
the solution at the exactly solvable Toulouse point. The 
quantitative agreement between their result and the TD- 
NRG for isotropic Kondo model is very good. However, 
our calculations also show that the transients towards the 
ultra long time behavior t/tK — * oo at time scale t ss tK 
are still dependent of the anisotropy ratio a = J z /J±. 
Deviations from a = 1 lead to deviations from an master 
curve for an isotropic Kondo model. Only at t — ► oo, we 
expect the merging of the curves. 
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FIG. 6: Spin-expectation value S z (t) for the isotropic Kondo 
spin with fixed J/D = 0.15 for different temperatures T = 
250, 33, 4.3, 0.075Tk and an initial magnetic field of H°/D = 
1 and J° = Jl = where T K /D = 0.00054. The lower panel 
(b) displays the same data as in (b) but plotted as function 
of t * T. NRG parameters: as in Fig[S] 




FIG. 7: Spin relaxation for the isotropic Kondo spin with 
fixed J/D = 0.15 for different strength of the final magnetic 
field Hi = 0.01, 0.1, 1.0, 10, 100T K for T = 1.2 * 10" 8 starting 
from H z /D = 0.1. The upper panel (a) shows the case of 
J° = 0, (b) the case J° = J 1 . The inset in (b) display an 
estimate for a spin-relaxation scale tspin as defined in the text. 
NRG parameters: as in Fig. 



2. Finite temperature relaxation 

Starting from a decoupled, fully spin-polarized impu- 
rity, the temperature dependence of the spin relaxation 
is plotted in Fig. for the isotropic, antiferromagnetic 
Kondo spin. The characteristic spin-relaxation time ti ow 
saturates for temperatures T < Tk and is of the order of 
1/Tk at very low temperatures, similar to the findings 
in the single impurity Anderson models The upraise of 
S z (t) for very long times at the temperatures T > Tk 
is an artifact of the insufficient time resolution at high 
temperatures. To illustrate this point, the same data is 
plotted as function of t * T in the lower panel (b) . 



3. Magnetic field dependence of the spin relaxation 

The time evolution of a polarized decoupled Kondo 
spin with Hamiltonian with a remaining finite local mag- 
netic field H 1 is plotted in Fig. {7\ For t < the lo- 
cal spin is fully polarized by a strong magnetic field. 
The figure shows the evolution of S z (t) for five differ- 
ent magnetic field strength H 1 = 0.01,0.1, 1, 10, IOOTr-, 
where the Kondo scale is given by Tk = 0.000534Z? for 
J/D = 0.15. For the upper panel (a), the local spin was 
decoupled from the bath for times t < 0, while the curve 
in the lower panel (b) were calculated for a Kondo spin 
coupled to the bath at all times. The difference in the 
time evolution between the two scenarios are very small 
in contrary to our previous study of the single impurity 
Anderson model 15 . The strong magnetic field destroys 
the Kondo correlation and the regimes are described by 
a spin polarized local moment fixed point in both cases: 



the Kondo spin-flip scattering is an irrelevant operator 
due to the external field. In the Anderson model, how- 
ever, the fixed points of both scenarios are completely 
different: in case of an absent hybridization, we have a 
free orbital or local moment fixed point while for finite 
coupling we started from a spin polarized mixed valent 
fixed point. 

Fig. clearly shows a decrease of the characteristic 
time scale with increasing final magnetic field. In order 
to quantify this behavior, we asked how long does it take 
until the spin-expectation value reaches some fraction of 
its asymptotic value S z (oo). For this purpose, we intro- 
duce the reduced spin expectation value function 



/(*) 



S z {t)-S z (oo) 
S z (0)-S z (oo) 



(73) 



and define a magnetic field dependent spin-relaxation 
time scale by the condition: f(t sp i n ) = 0.15. Note that 
the time scale t sp i n does not have the meaning of a re- 
ciprocal relaxation rate, since f(t) is - in general - not 
an exponential function. The choice of 0.15 generalizes 
the criterion given for ti ow at the end of section IVI A II 
While the value is arbitrary, of course, but setting to 
0.2 or 0.3 yields the same qualitative picture as depicted 
in the inset of Fig. EJd). For very weak magnetic fields, 
H < O.ITk, the spin-relaxation time t sp i n (H) remains 
almost field independent while for H > Tk the relax- 
ation time rapidly drops and reaches ti ow * D = 4.0 for 
(a) and t spm * D = 9.4 at H} = IOOTrt- 

In order to demonstrate the usefulness of the scale 
tspin(H), we rescale the data S z (t) of Fig. EJi) in dimen- 
sionless units using Eq. I|73|l and plot the resulting /(£) 
versus the dimensionless time scale t/t sp i„ . The result is 
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FIG. 8: Data of Fig.|7] panel (a) plotted as rescaled spin evo- 
lution /(£) defined in Eg. 1731 versus t 3p i„(H) for five different 
field strength and J° = 0. 



depicted in Fig. [S] For magnetic fields much larger than 
the Kondo temperature Tk , the spin relaxation is "fast" 
since the asymptotic value S z (oo,H) is reached rather 
rapidly. No Kondo correlations develop and no univer- 
sality is observed. For H\ < Tk, however, all curves 
collapse onto a universal curve f(t/t sp i n ) as clearly visi- 
ble in Fig. El 



B. The ferromagnetic regime 

1. Short and Long time Behavior 

As discussed in section IV IJI the ferromagnetic regime 
of the Kondo model is governed by the vanishing of J± for 
T — > 0. For the isotropic Kondo model, the Hamiltonian 
flows towards the local moment fixed point A2£ 

We start with a Kondo spin, always coupled to the 
bath, i.e. J a — constant and only subject to a switched 
magnetic held. In this case, the strong correlations al- 
ready form at times t < 0. In Fig. Et): S z {t) normalized 
to 52(0) is plotted for four different temperatures and two 
values of ferromagnetic couplings J = —0.1,-0.2. The 
spin-relaxation time increases with decreasing tempera- 
ture, since the coupling J(T) between conduction band 
and local spin is significantly reduced as predicted by the 
pour man scaling result l|65|) . In Fig. Et>, the same data 
as in panel (a) is plotted versus the dimensionless vari- 
able x = 2ppJ(T)Tt. All curves, which cover about six 
orders of magnitude of time dependency, collapse very 
well onto a single scaling curve. Only the effective value 
of J(T) entered the scaling variable x in addition to ex- 
plicit temperature. We did not need to include an time 
dependency in J(T) as suggested in Ref. |3^ 

The TD-NRG does fail to predicted the vanishing 
of spin-expectation value. We observe a saturation at 
approximately S z (0) * 0.55 which occurs at x ~ 1. 
Note, however, that x = 1 corresponds to t * T — 
[2pF J(T)] -1 3> 1 very far outside of the validity range of 
the method. The NRG does not have access to energies 
much lower than the temperature T. Nevertheless, the 
scaling predictions of the TD-NRG gives very strong in- 
dications that the spin-relaxation rate is proportional T 
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FIG. 9: Normalized spin-expectation value for the isotropic 
Kondo model in the FM regime for two different values 
J/D = —0.2, , —0.1, and four different temperatures T/D = 
0.13,0.018,0.0023,4* 10~ 5 . While J = J z = J ± is kept 
constant, a magnetic field of H z /D = 0.01 is switch off at 
t = 0. Fig (a) shows the time evolution of the normal- 
ized spin expectation value S z (t) = (S z )(t) while in (b) the 
same data is plotted versus the dimensionless scaling vari- 
able x = 2pfJ(T)T * t. The result is correct only for x < 1 
since for large x energy scales smaller than the smallest energy 
scale accessible for the NRG at a give temperature T would 
be needed^. Parameters: A = 1.5; N 3 - 1000, a d - 0.1, 
N z = 16. 



for the isotropic ferromagnetic Kondo model in addition 
to the expected dependency of the reduced Kondo cou- 
pling J(T). At low temperature, the effective coupling 
renormalized to zero and, therefore, the perturbative ap- 
proach of Langreth and Wilkins* gives the correct an- 
swer for the relaxation rate T ex (pFJ) 2 krjT in the long 
time limit. The proportionality of F and the tempera- 
ture stems from the fact that the number of hnal states 
for the scattered electrons is proportional to the thermal 
spread fc^T. 



2. Spin precession 

So far, we have only focused on the spin relaxation from 
a spin-polarized state along the z-axis. Precession of the 
spin around an external held requires a magnetization 
perpendicular to this external held. To accomplish this, 
we applied a hnite magnetic held H = Q.lDe x for t < 
to induce a magnetization in ^-direction. At t = 0, the 
external held is rotated in z-direction. In the absence of 
any relaxation mechanism, i.e. for J = 0, the spin would 
only precess in the xy-plane once H is aligned along the pr- 
axis. The hnite Kondo interaction, however, induces spin 
flips and causes spin exchange between the conduction 
band and the localized spin. In the thermodynamic limit 
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FIG. 11: Relaxation rate for the S^ and the S z component 
of the spin by fitting S x (t) — 0.5cos(u>Lt) exp(— T x i) and 
S z (t) = S z (oo)[l — cos(u;£t)exp( — T z i)] for isotropic ferro- 
magnetic and anti-ferromagnetic coupling. NRG Parameters: 
N s = 1000, iV 2 = 16, A = 1.8, Niter = 60, a d = 0.1. 



FIG. 10: Spin precession of a Kondo spin coupled to the con- 
duction electrons at all time J z ~ J± = const for two dif- 
ferent values of J. At t — the magnetic field H x /D — 0.1 
is rotated around the t/-axis and pointing along the z-axis, 
H = 0.1De z for t > 0. NRG Parameters: N s = 1000, N z = 
16, A = 1.6, Niter =60, cm = 0.1. 



and for t — ► oo, the expectation value of the spin vector 

(5) (t) will evolve into a spin aligned along the z-axis. 

The time evolution of all spin components is displayed 
in Fig.^]|for a magnetic field of constant strength rotated 
from the x- to the z-axis at t = 0. Since the Zeeman term 
of the Hamiltonian (|59|l does not commute with the to- 
tal S z and only with Sf ot , the TD-NRG calculations are 
much more time consuming. The black and red curve 
(online) show the S x (t) and S y (t) components, respec- 
tively. As expected, the spin precesses in the xj/-plane 
which is clearly visible for J/D = 0.2. In addition, these 
two components decay with time and vanish for long 
times, while the S z (t) component grows and reaches a 
value close to the thermodynamic expectation value with 
respect to Ti* . This indicates clearly that the TD-NRG 
can covers more complex relaxation phenomena. Indeed 
the local expectation values tent to approach their long 
time thermodynamic limit, which is a highly non-trivial 
feature of the method. 

For a quantitative analysis, we fitted the two spin com- 
ponents S x (t) and S z (t) by the following phcnomeno- 
logical ansatz: S x (t) — S x (0) cos(ujl) exp(— T x t) and 
S z (t) = 5 z (oo)[l — cos(u> z t) exp(—r z t)]. S z (t) shows an 
oscillatory behavior with a frequency lo z l rs 0.006 for the 
ferromagnetic Kondo coupling and an uj z l rs 0.02 for the 
anti-ferromagnetic Kondo coupling. The decay of the 
local S x component leads to a very small but finite mag- 
netization in the conduction electron chain in x direc- 
tion. Its local component induces a spin precession of 
the S z which is damped by the spin-flip processes. For 
a chain length of N = 60, we estimate the average mag- 
netization per chain link as S x (0)/N = 1/120 w 0.0083 
which is a good estimate for the lower limit of an ef- 



fective local magnetic field in x direction generated by 
the conduction electrons in a finite size system. In the 
anti-ferromagnetic case, the coupling of the local spin 
to the conduction electrons is strongly enhanced which 
might cause an inhomogeneous distribution of magneti- 
zation over the Wilson chain. This effect must vanish 
in the thermodynamic limit N — ► oo,A — > 1 + . We find 
(T X ,T Z ) = (0.014,0.019) for J/D = 0.2 and (0.12,0.113) 
for J/D = 0.4. 

In Fig. ^2 the effective spin-relaxation rates T x and T z 
are displayed as function of \J/D\ for an initially decou- 
pled local spin, i.e. J(t < 0) = 0. With the exception of 
J/D = 0.25, T x is smaller than T z for anti-ferromagnetic 
coupling while the transversal relaxation is always slower 
than the longitudinal for ferromagnetic coupling. How- 
ever, the absolute values strongly depend on the fitting 
functions for S x (t) and S z (t), so that the relaxation rates 
are roughly equal for all directions. Note that the usual 
dephasing time T2 in NMR and ESR experiments origi- 
nates from an average over signals stemming from many 
spins subject to local fluctuations of the externally ap- 
plied field. This inhomogeneities yields the ESR condi- 
tion I/T2 > 1/Ti which does not necessarily hold for the 
single spin decay investigated here. 



VII. SUMMARY AND DISCUSSION 

A. Discussion 

We have presented a powerful novel approach to the 
nonequilibrium dynamics of quantum impurity systems 
based upon Wilson's NRG. Its virtue is based on three 
key points: (i) all states of the Wilson chain are used for 
the time-evaluation of local operator, (ii) time expecta- 
tion values can be expressed by a summation of all NRG 
iterations, where at each iteration the matrix elements 
of any local operator is weight by the reduced density 
matrix which contain information about dissipation and 
decoherence and (iii) the continuum limit of the bath is 
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mimicked by using the z-trick where the true time evo- 
lution is obtained by averaging over a different bath dis- 
cretization. By its nature, it is applicable to arbitrary 
temperatures, and allows for the study of the tempera- 
ture evolution of real-time dynamics. 

We have previously established the accuracy of the ap- 
proach on all time scales up to 0(tT) by comparison 
with exact solution of the RLM 15 . In this paper, we 
additionally benchmarked our method to the exact an- 
alytical solution of the decay of the off-diagonal matrix 
element poi of the local density matrix which is taken 
as measure for decoherence of a spin state. This is the 
first application of the TD-NRG to bosonic baths. In 
this case, the TD-NRG is much less susceptible to dis- 
cretization errors than for fermionic baths. This behavior 
roots in the extended number of bath states, in principle 
infinite, which are available for distributing energy and 
phase at each energy scale D^. From this observation, we 
conjecture that the TD-NRG will be even more suitable 
from bosonic baths than for fermionic baths. A detailed 
study of the real-time dynamics of the spin-boson model 
in the subohmic regime will be published soon. Com- 
bining fermionic and bosonic bathsai will open new pos- 
sibilities for the description of spin decay of ultra-small 
quantum dots coupled to the leads as well as charge noise. 

We have investigated the spin decay of the Kondo 
model in the ferromagnetic and in the antiferromagnetic 
regime. For the ferromagnetic regime, we identified a 
universal dimcnsionless time scale x — 2ppJ(T)Tt for 
the transient behavior of the spin decay. In the antifer- 
romagnetic regime, the spin response at short time scales 
follows excellently the presented analytical perturbative 
solution which is valid up the second order in the coupling 
constant. At long times, the Kondo correlations domi- 
nate, and the rescaled curves appear to collapse onto a 
universal curve for t/tx — > oo as predicted by conformal 
field theory^. At short and intermediate times t ~ tj<, 
the anisotropy of J z and Jj_ influences the time evolution 
and no universality is found as a one parameter scaling 
for different values of J z /J±. This reflects the different 
role of J z and J± in the spin-relaxation process. The 
spin relaxation is governed by the spin-flip term Jj_ and 
not by J z which only enters through renormalization in 
higher order processes into J± . Only at T — and at in- 
finitely long times, the anisotropy Kondo model becomes 
isotropic in the strong coupling fixed point. Here, we 
expect the validity of universal scaling. 



B. Outlook 

The time-dependent NRG approach will open new 
doors to our understanding of a new class of nonequi- 
librium problems whose real-time dynamics is governed 
by a strong entanglement of the environment with the 
quantum-impurity states. The exact analytical solution 
of decoherence of a pure quantum states explicitly shows 
that the true dynamics is not simply given by an expo- 



nential decay rate T but a function T(t). The reproduc- 
tion of these non-trivial analytical results of the spin- 
boson model with all details reveals already the strong 
potential of our method. The application to charge- 
transfer reactions in biological systems 2 - influenced by 
the molecular vibrations is subject of an ongoing research 
project. We hope that this will ultimately lead to a bet- 
ter understanding how secondary and tertiary molecular 
structures of protein molecules influence the chemistry 
and, therefore, the functionality of such proteins. It was 
pointed out by Schulten^i that the chemical structure of 
reaction centers in protein molecules by itself does not 
give us the desired understanding why thermal fluctua- 
tions do no harm the deterministic outcome of the com- 
plex cellular chemistry which is the basis of existence of 
living organisms. 
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APPENDIX A: CALCULATION OF THE 
OVERLAP MATRIX ELEMENTS 

In this appendix, we describe in detail the calculation 
of the overlap matrix elements of Eq. (|34fl . We need the 
matrix element 



S rtS (m) = (s;m|r;m) 



(Al) 



in order to calculate the reduced density matrix in basis 
of the Hamiltonian Ti^ n needed in Eq. i|35|) . Here, \s;m) 
is an eigenstate oiTL l m and \r; m) an eigenstate of Ti^ and 
the environment variable e has been droped for simplicity. 

i/f 

Independent of the dynamics of Tim , the original basis 
of the Hamiltonian matrices H* and EL is identical prior 
to diagonalization. Therefore, the overlap is given by 






(a imp ; 0|a- 0) = 5 at a ' ( A2 ) 



where ai mp labels the states on the impurity. Let us 
assume, we know all matrix elements S r ^ s (m) at iteration 
to. Then, the NRG recursion relation adds the same 
degrees of freedom a m+ \ to the new chain of length to, 
independent of the dynamics of "H„ l+ i- Since the new 
eigenstates can be expanded into 



\r:m + l) = ^ p lki a in+i]\k, "m+i) (A3) 

|s;m + l) = Y^ Pik'[ a m+i}\k',a m+1 ) , (A4) 



Otm+l,k 
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the matrix element S r '_ s '(m + 1) is given by 
S r >,„>(m + 1) = (s';m + l|r';m+l) (A5) 

k,k' a rri j r i .a' , - 

x{k,a m +i\k',a' m+1 ) 

k,k' a m+1 



xp:> k > 



}Sk,k'(m) 



In combination with the initial condition (|A2|) . Eq. (|A6|) 
defines a recursion relation from which all overlap matrix 
elements are obtained. 
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